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Abstract 


A systematic study is presented to guide the selection of a numerical solution strategy for 
URANS computation of a subsonic transport configuration undergoing simulated forced 
oscillation about its pitch axis. Forced oscillation is central to the prevalent wind tunnel 
methodology for quantifying aircraft dynamic stability derivatives from force and moment 
coefficients, which is the ultimate goal for the computational simulations. Extensive 
computations are performed that lead in key insights of the critical numerical parameters 
affecting solution convergence. A preliminary linear harmonic analysis is included to 
demonstrate the potential of extracting dynamic stability derivatives from computational 
solutions. 


Nomenclature 


C L 

= 

lift coefficient 

C m or CM 

= 

pitching moment coefficient 

Cref 

= 

mean aerodynamic reference chord, =10.98 in. 

f 

= 

sinusoidal oscillation frequency about pitch, roll, or yaw axis, Hz. 

k 

= 

reduced frequency, =2jt f • (c re /2)/U CX) 

log(r/r0) 

= 

order of magnitude drop in solution residual, e.g. -3 means 3 orders magnitude reduction 

A4 US M3 D 

= 

USM3D freestream Mach number, =0.2 

Moo 

= 

experimental freestream Mach number, =0.077 

n, n+1 

= 

indices for numerical time step 

R^cref 

= 

Reynolds number based on c re f 

ratel 

= 

primary VGRID “viscous” stretching factor, see Eq. 1 

rate2 

= 

secondary VGRID “viscous” stretching factor, see Eq. 1 

At 

= 

incremental time step, sec 

. * 

At 

= 

characteristic time step, =At 0 U oo /c re f 

U CO 

= 

freestream velocity, =86 ft/sec 

H" 

y 

= 

turbulence wall coordinate 

a 

= 

angle of attack (static: a, deg, dynamic: a=a 0 +Aa, deg) 

a 0 

= 

nominal angle of attack during pitch oscillation, deg. 

Aa 

= 

incremental angle of pitch oscillation about body axis, deg. 

& i+i 

= 

VGRID “viscous” grid spacing normal to surface at node /+ 7, see Eq. 1, inches 

<5; 

= 

spacing of first node off of surface (i=l) in “viscous” grid layers, see Eq. 1, inches 
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Key Acronyms 

GTM 

SA 

SST 

URANS 

14X22 


Generic Transport Model 

Spalart-Allmaras one equation turbulence model 

Menter’s Shear-Stress Transport two equation turbulence model 

Unsteady Reynolds-Averaged Navier-Stokes 

NASA Langley 14- by 22-Foot Low-Speed Wind Tunnel 
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Figure 1. Fatal accident distribution for 
commercial transports and general aviation. 
Source NSTB database 1990-96 


I. Introduction 

The NASA Aviation Safety Program is sponsoring a multi-year effort through the Integrated Resilient Aircraft 
Control (IRAC) project to improve the survivability of civil transport aircraft that encounter critical loss-of-control 
(LOC) flight upset conditions. IRAC’s intended goal is the creation of validated, multidisciplinary integrated aircraft 
control design tools and techniques that can be used to enable safe flight in the midst of adverse flight conditions 
(i.e. faults, damage, icing, and/or upset) by addressing aircraft stability and control (S&C) parameters [1]. 

In the late 1990’s the US Government along with Industry formulated a partnership to address and develop plans 
to reduce aviation accidents and loss of life. That partnership developed into the Commercial Aviation Safety Team 
(CAST). During their investigation of the data it was learned that the highest contributor to loss of life were LOC 
events as can been seen in Fig. 1 from data gathered by the National Safety Transportation Board (NSTB) [2]. 

A series of intervention strategies were developed to 
reduce the occurrence of LOC accidents, which can arise 
from a variety of causes. From the initial analysis and 
developmental phase it was determined that the current 
simulation technology was inadequate to train and prepare 
pilots for realistic conditions outside of the normal flight 
envelope (i.e. stall, post-stall and upset conditions). 

Traditionally there have been three methods of choice used in 
analyzing S&C characteristics. Flight-testing the actual 
aircraft, would be the obvious first and most accurate, though 
not necessarily, most economical choice [3,4]. This method is 
prohibitive not only because of the cost involved, but also 
because of the inherent danger posed to the flight crew as a 

result of an inadequate escape mechanism in transport aircraft for cockpit members. Wind tunnel testing on scaled 
models is the second method of choice, but it too has large cost and time requirements [3]. In addition to cost and 
time concerns, other factors such as Reynolds-number scaling, blockages and interference to the flow caused by the 
model support system prevent proper modeling of the full scale behavior and contribute to limiting the use of this 
technique. Incorporating the use of data sheets, linear aerodynamic theory, and empirical relations in the third 
method can be effective in typical flight maneuvers but is very limited in its ability to predict flight dynamics at the 
outer edges of the flight envelope in stall, post-stall or upset conditions [3]. 

Figure 2 portrays a graphical representation of the 
limitations of the current state of the art in predicting flight 
dynamics in the outer edges of the flight envelope. The 
chaotic red line, representing “black box” data from an 
aircraft experiencing a LOC accident, reveals a flight 
trajectory well outside of the available wind tunnel data. Two 
of the main reasons behind the limited flight data are: 1) the 
current data sets meet the minimum requirements for 
simulator certification, which do not typically explore beyond 
normal flight conditions and 2) flight and wind-tunnel test do 
not usually provide aerodynamic measurements of upset 
conditions. To bridge this data gap, Computational Fluid 
Dynamics (CFD) is being utilized to analyze flight dynamics 
in the outer reaches of the flight envelope. 

This study contributes to the Integrated Dynamics & 

Flight Controls (IDFC) subproject within IRAC. Within the modeling components IDFC, the goal is to develop 
reliable computational tools for predicting and analyzing stability and control characteristics of aircraft in damaged 
and upset flight conditions. A prior study reported in Ref. [5] advances this objective by assessing the utility of two 
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Figure 2. Graphical 
available aerodynamic data for flight envelope 
undergoing LOC [2]. 
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unstructured Navier- Stokes flow solvers for capturing the degradation in static stability and aerodynamic 
performance of a subsonic transport due to airframe damage. The present work delves into the complex realm of 
dynamic stability in upset flight conditions. The focus is to assess the applicability of CFD modeling for dynamic 
aircraft motions on the post-stall aerodynamics that are characteristic of LOC flight. The following study has been 
initiated to develop the practices for computing the aerodynamic and stability characteristics of an aircraft 
undergoing periodic oscillation about a post-stall angle of attack using the tetrahedral-based USM3D unsteady 
Reynolds Averaged Navier-Stokes (URANS) flow solver. To be assessed are 1) the accuracy and efficiency of the 
NASA USM3D flow solver for predicting the hysteresis behavior of force and moment coefficients for a subsonic 
transport model undergoing pitch oscillation, and 2) the sensitivities of USM3D input parameters affecting time- 
accuracy of the aircraft in dynamic motion which will contribute to the development of practical guidelines for 
applying USM3D for future dynamic S&C problems. This effort supplements and reinforces the findings from a 
similar study performed on unmanned military air vehicle [6]. While it is recognized that URANS methodology may 
not be suitable for massively separated flows, the intent is to lay a solid numerical foundation for future studies, 
which will include more advanced turbulence models. 

II. Geometry and Experiment 

The IRAC focus configuration (Fig. 3) is a representative twin-engine commercial transport configuration called 
the Generic Transport Model (GTM) that has received extensive static and dynamic force and moment testing in the 
NASA Langley 14- by 22-Foot Subsonic Tunnel (14X22) [7-9]. Conventional forced-oscillation tests conducted in 
2001 and 2003 were expanded upon in 2007 and 2008 to include a broader range of oscillation frequencies that 
would allow identification of unsteady aerodynamic mathematical models. The model span is 82.2 inches, and the 
mean aerodynamic reference chord, c re /, is 10.98 inches. The tests were conducted at sea level atmospheric 
conditions with a dynamic pressure of 10 lb/ft 2 , corresponding to a freestream velocity U =86 ft/sec, a Mach 
number Moo =0.077 and chord Reynolds number Re cre f= 0.54xl0 6 . Boundary layer transition strips of width 0.1 inch 
and number 36 grit were applied at 10% chord on the upper and lower surfaces of the wing, horizontal and vertical 
tails. A 0.1-inch wide ring of number 36-grit was applied 1 inch aft of the apex of the nose (measured 
longitudinally). This gritting pattern was chosen to maintain commonality with previously tested configurations. At 
the time of this writing it is unclear if this technique was effective in tripping the flow (there was no verification 
performed, and with the exception of the nose grit, this particular tripping technique differs from well established 
low-speed practices, e.g. Braslow [10]). However, these are the limitations of dynamic wind tunnel testing that 
might be expanded by CFD application. 

Forced oscillation tests were conducted for the baseline and damaged GTM geometries. The configurations were 
oscillated over a nominal range of roll, pitch, and yaw body angles in a sinusoidal fashion at frequencies and 
amplitudes that enable the identification of both conventional aerodynamic characteristics and advanced 
mathematical models for rigid-body aerodynamics in nonlinear unsteady flight regimes. Oscillation tests were 
conducted individually in each axis to acquire aerodynamic damping characteristics. For this study sinusoidal pitch 
oscillation of Aa =± 5 deg was implemented about a 0 = 24 deg with frequencies /= 0.43, 0.12, and 0.05 Hz, which 
correspond to reduced frequencies k=2jzf-( c re f/2)/U =0.01438, 0.00401, and 0.00167, respectively. 

III. Computational Method 

The cell-centered computations are produced from the USM3D Navier-Stokes flow solver [11] that is part of the 
NASA Tetrahedral Unstructured Software System (TetrUSS) [12]. USM3D is a parallelized tetrahedral cell- 
centered, finite volume Navier-Stokes flow solver. The term “cell centered” means that the finite volume flow 
solution is solved at the centroid of each tetrahedral cell. Inviscid flux quantities are computed across each 
tetrahedral cell face using various upwind schemes. Spatial discretization is accomplished by a novel reconstruction 
process, based on an analytical formulation for computing solution gradients within tetrahedral cells. The solution is 
advanced in time [13] by a 2 nd -order “physical” time step scheme, a 2 nd -order “dual” time step scheme, or to a 
steady-state condition by an implicit backward-Euler scheme. Several turbulence models are available: the Spalart- 
Allmaras (SA) one-equation model, the two-equation Jones and Launder ks turbulence model, the Menter Shear 
Stress Transport (SST) two-equation model, and the nonlinear Algebraic Reynolds Stress Models of Girimaji and 
Shih/Zhu/Lumley. Detached Eddy Simulation has been implemented in all of the turbulence models. A capability to 
trip the flow at specified locations on aerodynamic surfaces has been implemented for the k - 8 turbulence model. 
USM3D has capabilities for dynamic grid motion, which is being utilized in the current study. 
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IV. Grid Generation 

Half-span tetrahedral grids of nominally 6-, 12-, and 24- million cells were generated using the unstructured 
tetrahedral grid generation codes GridTool and VGRID, which are part of the NASA TetrUSS system. The surface 
triangulation for the 6-million cell grid is shown in Fig. 4. A developmental version of VGRID [14] was used for the 
present study. Thin-layer tetrahedral grids were generated to meet requirements for cell-centered computations from 
the USM3D flow solver. The near-wall spacing prescribed to achieve a turbulent wall coordinate y + of 0.75 in the 
tetrahedral cell center near the mid-chord of the c re f. Since the VGRID Advancing Layers Method [15] marches 
nodes away from the surface (which are subsequently connected to form tetrahedral cells), then an initial VGRID 
“viscous” spacing corresponding to a y + =3 of the first node was prescribed in order to achieve the y + =0.75 at the 
first cell centroid. For the GTM at the wind tunnel test conditions, the required initial first-node wall spacing is 
8 1=0. 001239 inches, and stretching factors of ratel=0.15 and rate2=0.02 , where the nodal spacing layers are 
defined by the Eq. (1). 


d i+1 =d 1 [\+ratel(\+rate2) 1 ] 1 (1) 

This spacing distribution results in approximately 54 tetrahedral cells (18 nodes) across the boundary layer at the 
mid-chord of the c re f. The triangular surface faces were stretched as high as a 10-to-l ratio along the leading edges of 
the wing and tails. 


V. Characterization of GTM Aerodynamics 

The character of the static lift and pitching moment polar is presented in Fig. 5. The 14X22 static wind tunnel 
data [7] extends up to 85 deg angle of attack. Computed results are included at angles of attack of 4, 10, 14, 20, 24, 
35, and 40 degrees on half-span grids of 6-, 12-, and 24-million cells. The static USM3D solutions are computed as 
time-accurate URANS flow solutions assuming fully turbulent flow with the SA and SST turbulence models at a 
characteristic time step of At*=At 0 U oo /c re f =0.02 for 1500 time steps. An assumption is made that the wind tunnel 
sting and walls can be ignored in these calculations. 

Sensitivity to grid density is portrayed in Fig. 5 for half-span grids of 6-, 12-, and 24-million cells using the SA 
turbulence model. Both the lift and pitching moment coefficients in Fig. 5 indicate a break in slope around a =10 
deg that is indicative of a wing-stall onset. Lift continues to increase up to a maximum around a=35 to 40 deg, 
beyond which the aircraft is in deep stall. The pitching moment continues to lose pitch up authority with increasing 
angle of attack. The computational solutions suggest that the effect of grid density is minimal relative to their 
correlations with the wind tunnel data with the exception of o=14 deg, where the 6-million cell grid encounters 
significant differences relative to the 12- and 24-million cell grids. The computations generally over-predict the lift 
and under-predict the pitching moment in the post-stall region, which is an expected result due to excessive levels of 
turbulent viscosity produced at those conditions by the SA turbulence model. 

The impact of turbulence model is addressed in Fig. 6 for the 6-million cell grid. The dominant effect is in the 
stall region beyond a=10 deg. Predicted lift for the SST model is significantly below that of the SA model and 
experimental data over angle of attack range from 14 to 24 degrees. Significant differences persist at the higher 
angles of attack with the SST model correlating best at cc= 40 deg. The SST model is known to produce lower levels 
of turbulent viscosity than the SA model for massively separated thereby allowing disturbances to propagate more 
readily to the wing upper surface. 

The focus of the present study is on establishing guidelines for computing dynamic S&C characteristics by 
mimicking the experimental approach of imposing a periodic oscillation about a nominal angle of attack. Results 
from the dynamic 14X22 wind tunnel data [9] are presented in Fig. 7 for a 0 = 24 deg with Aa=+/-5 deg of pitch 
oscillation. The pitching moment coefficients are plotted relative the oscillating body angle, Aa. Both raw and mean 
(averaged) experimental data are presented over 40 pitch cycles for oscillation frequencies of 0.43, 0.12, and 0.05 
Hz. The scatter of the raw data (plotted in light gray) reflects the difficult nature of experimental dynamic testing 
where both the aerodynamic behaviors of interest as well as systematic and random measurement errors are present. 
This data is averaged to create the “mean” wind tunnel data that will be used for comparisons in the computational 
study. 
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VI. Harmonic Analysis and Model Identification 

Conventional flight dynamics simulations incorporate the aerodynamic model by assuming the aerodynamic 
forces and moments can be approximated by a series expansion in terms of aircraft states and controls. The terms in 
the expansion are the stability and control derivatives. A large database of knowledge has been built around these 
derivatives providing substantial insights for analysis and design. Several experimental methods can be used to 
obtain derivatives with respect to rates (damping) but the conventional method is to use forced-oscillation tests. 
Conventional analysis of longitudinal forced-oscillation data assumes that the aerodynamic coefficient is a linear 
function of angle-of-attack, pitching velocity, and their rates. This approach [16] produces in-phase and out-of-phase 
coefficients, C Qa and C a ^ respectively, where subscript “ a ” represents either normal force, C N , or pitching moment, 

C m . A method of harmonic analysis [17] can be applied to the measured aerodynamic forces and moments to allow 
estimation of the in-phase and out-of-phase coefficients. For harmonic analysis a mathematical model is postulated 
as 

r r 

C a (t) = A 0 + ^ Aj cos (jcot) + ^ Bj sin (jcot) a =N or m (2) 

7=1 7=1 

where A 0 , Aj , and Bj are the Fourier coefficients. For a model with linear aerodynamics, the in-phase and out-of- 
phase components of C a (t) can be expressed in terms of the 1 st harmonic Fourier coefficients as 



= — = C ~k 2 C 

a =N or m 

( 3 ) 


n a <? 

a a 





A, _ _ 



c a 

= -^ = C a +C « 

a =N or m 

( 4 ) 

U q 

ka A q 



Also shown in Eqs. (3-4), are the commonly expressed relationships of these coefficients to the steady flow 
damping terms and acceleration terms. When an aircraft is operating at low angles-of-attack or in linear and steady 
aerodynamic flight regimes the out-of-phase component is a good approximation for the damping. This formulation 
works well with the series expansions used in flight simulations. However, at high angles-of-attack or regions where 
separated flows dominate, the acceleration terms can become important and demonstrate strong frequency 
dependence. Under these conditions the nonlinear unsteady aerodynamic behaviors are not well captured by this 
term alone and more advanced models are required. 

Although the conventional model structure has limitations, substantial insights can be gained by comparing the 
in-phase and out-of-phase components derived from wind tunnel tests against that obtained from CFD calculations. 
The first harmonic is an important term even under conditions where higher order terms might be required. Due to 
the orthogonal nature of the model terms the first harmonic is not changed even when higher order terms are 
required to model nonlinear aerodynamic behavior. Also, the first harmonic will demonstrate frequency dependence 
in the presence of unsteady aerodynamics. Thus two key aspects of the aerodynamic response can be captured by 
comparing the aerodynamic force and moment time histories and the 1 st harmonic coefficients between wind tunnel 
and CFD results: (1) time history comparisons demonstrate that full nonlinear aerodynamic responses have been 
captured, and (2) the 1 st harmonic coefficients demonstrate the levels of damping and unsteady aerodynamics if 
present. 

VII. Approach to Numerical Guideline Study 

The primary goal for this paper is to establish guidelines for applying USM3D toward computing reasonably 
accurate dynamic force and moment characteristics on the GTM aircraft with as little computer resource as possible. 
An approach similar to that outlined in Ref. [6] is applied to explore the relationship between the numbers of time 
steps and subiterations of the 2 nd -order time step scheme. The bulk of this work will focus on understanding the 
relationships between numbers of time steps per oscillation cycle and subiterative convergence between each time 
step to establish application guidelines for computing dynamic stability characteristics for a spectrum of cases. 

The prevalent experimental methodology for analysis of dynamic stability characteristics of an aircraft calls for 
extracting the stability derivatives from the force and moment coefficients as the aircraft undergoes a sinusoidal 
motion about an axis, as illustrated in left image of Fig. 8. A similar approach is mimicked in the present work for 
longitudinal pitch oscillation. 

The strategy from [6] incorporates visual correlation to understand the relationship between the number of time 
steps per oscillation cycle and subiterative convergence at a time step representative of difficult convergence. As 
illustrated in Fig. 8, a set position within the pitch oscillation cycle (in this case at -3 deg downswing) is selected to 
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monitor the correlation between the pitching moment coefficient and its subiterative convergence between the 
accompanying time step n and n+1. The middle plot in Fig. 8 depicts a dynamic pitching moment coefficient, C m , 
that exhibits hysteresis as the GTM oscillates +/- 5 deg about its pitch axis at a 0 =24 degrees. (The black dashed line 
is “mean” wind tunnel data from [9].) The arrow leaning to the right denotes the upswing, and the arrow leaning to 
the left denotes the downswing. The right-hand plot in Fig. 8 conveys the corresponding order-of-magnitude 
convergence in the USM3D solution residual drop ( log(r/rO )) highlighted in red, and C m highlighted in blue. In this 
example, the solution residual has dropped over 3 orders of magnitude and the C m is converged after 100 
subiterations using 360 time steps per pitch oscillation cycle 

This study is organized by first establishing a set of “well converged” reference solutions. These solutions will 
be generated by applying 2 nd order “physical” time stepping in USM3D with the Spalart-Allmaras (USM3D/SA) 
turbulence model in a manner that balances the numbers of time steps per pitch oscillation with the number of 
subiterations between time steps, to maintain a fixed number of total iterations. A sufficiently large number of 
subiterations will be used to insure a well-converged solution between time steps. Next, reducing the numbers of 
subiterations for a selected oscillation frequency and monitoring its subiterative convergence and resulting departure 
from the reference solution will identify an efficient solution strategy. Once a successful strategy is determined, it 
will then be tested with other pitch frequencies, the SST turbulence model, a larger pitch amplitude of Aa=+/-10 
deg, and for grid sensitivities. 


VIII. Results and Discussion 

The selected GTM test condition is a nominal post-stall a 0 = 24 deg with a pitch oscillation of Aa =+/- 5 deg at 
frequencies of ^0.43, 0.12, and 0.05 Hz that correspond to reduced frequencies &=0. 01438, 0.00401, and 0.00167, 
respectively. The flow conditions are M US m3d = 0.2 and Re cre f =0.54x1 0 6 . The USM3D results were generated at a 
slightly higher Mach number than the wind tunnel freestream 717^=0.077 to improve the numerical convergence with 
the assumption that compressibility effects are negligible. Fully turbulent flow is assumed for all solutions. The flow 
computations were performed with Roe’s Flux Difference Splitting [18] and no flux limiter. 

A. Establishing “Well Converged” Reference Solutions 

Table 1 shows the matrix of computational runs used to 
generate “well converged” USM3D solutions using the SA 
turbulence model at pitch frequencies of 0.43, 0.12, and 0.05 Hz. 

Time steps were chosen to yield 360, 720, and 1440 steps per 
pitch cycle. Corresponding levels of subiteration between time 
steps, 100, 50, and 25, respectively, were defined to yield a total 
of 36,000 solution iterations over each pitch cycle. The cells 
highlighted in blue will be used in the next section in search of 
more efficient solution strategy. 

Figure 9 depicts plots of sub-iterative convergence of solution 
residuals (red) and pitching moment coefficient (CM) (blue) 
between time steps n and n+1 near the downswing Aa= -3.5° 

(a=20.5°) for the cases in Table 1. The order-of-magnitude drop 
of the solution residual is denoted on the left axis, and CM is on 
the right axis. The axis increments of CM are held constant 
between plots. The most prominent trend in Fig. 9 is that as pitch 
frequency decreases (left to right), the order of magnitude 
residual drop increases from approximately -4 to the about -2.5 suggesting a less robust convergence, and the 
converged CM becomes more unsteady. At each frequency, however, the quality of the well-converged subiterative 
convergence is similar (top to bottom) regardless of the per-cycle [time step]/[subiteration] combination 360/100, 
720/50, or 1440/25. 

The dynamic pitching moment coefficients for the cases in Table 1 and Fig. 9 are plotted in Fig 10. The solutions 
exhibit almost no sensitivity of the C m hysteresis to these prescribed combinations of time step and sub -iteration, 
indicating that solution accuracy is predominantly a function of total iteration count. The general shapes of the C m 
hysteresis curves are remarkably similar between CFD and experiment. Recall from the discussion of Fig. 7 that the 
wind tunnel data was derived as an average or “mean” of highly unsteady dynamic measurements. This “mean” is 
similar to averaging that occurs due to the integration used with harmonic analysis methods to determine the in- 
phase and out-of-phase dynamic stability derivatives [16]. The CFD is computed with the assumption of “Reynolds- 


Table 1 - GTM cases for “well converged” 
solutions with constant 36,000 total 
iterations. USM3D/SA, /= 0.43, 0.12, 

0.05Hz, a 0 =24 deg, Aa =+/- 5 deg. 




Pitch freq, Hz 
(reduced frequency, k) 

Time- 
step / 
cycle 

Number 
Sub iter 

0.43 

(0.0144) 

0.12 

(0.0040) 

0.05 

(0.0017) 

360 

100 

360/100 

360/100 

360/100 

720 

50 

720/50 

720/50 

720/50 

1440 

25 

1440/25 

1440/25 

1440/25 
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Averaged Navier-Stokes” [19], which completely filters out any unsteady turbulent eddies in this massively 
separated post-stall flow. Thus, it is intriguing that these two fundamentally different approaches to averaging 
produce such similar results. Before pursuing more costly turbulence models that resolve turbulent eddies, such as 
Detached Eddy Simulation, further study of the URANS formulation might be warranted to explore its applicability 
and limitations for modeling post-stall dynamic stability characteristics. 

B. Searching for Efficient Solution Strategy 

Table 2 conveys the y=0.43 Hz cases where sub iterations 
are reduced consistently for each 360, 720, and 1440 time step 
strategy to maintain comparable levels of total iterations. The 
goal is to determine the least number of total iterations needed 
to maintain reasonably converged dynamic C m . The blue 
highlighted cells from Table 1 are repeated in Table 2. 

The convergence plots in Figure 1 1 are taken from the “well 
converged” solutions for J=0A3 Hz in Figure 9 and 10. Arrows 
at their respective positions on the convergence plots denote the 
reduced sub iterations from Table 2. The corresponding C m 
hysteresis loops are included in Fig. 11. From this, a visual 
correlation can be established between the subiterative 
convergence and the convergence of the C m hysteresis loops 
relative to the “well converged” solution. It can be observed in 
Fig. 11(a) that the shape of the C m hysteresis changes from the “well converged” solution more rapidly below 30 
subiterations with 360 steps per cycle. Hence reasonably converged dynamics solutions might be achieved on the 
GTM with as few as 10,800 total iterations per pitch cycle, which will be further tested for two lower frequencies in 
the next section using the highlighted cases (360/100, 360/50, and 360/30) in Table 2. 

The question of “how much accuracy is needed” remains unanswered and may be resolved using harmonic 
analysis and system identification techniques. These techniques add the requirement that sufficient sample rates be 
used to satisfy the Shannon Sampling Theorem. It is entirely possible that the accuracy with the 10 subiterations in 
Fig. 11(a) having 3600 total iterations per cycle is more than adequate for determining suitable dynamic stability 
derivatives for flight simulator models. The present study is primarily focused on establishing guidelines for 
numerical convergence. 

C. Testing Strategy at Varying Frequencies 

The next phase in the analysis involves testing 
the new strategy at lower oscillation frequencies 
of 0.12, and 0.05 Hz to assess its validity at other 
dynamic conditions. The highlighted cells for 
y=0.43 Hz from Table 2 are repeated in Table 3. 

Results for the cases in Table 3 using less 
subiteration are plotted in Fig. 12 where small 
differences are noted in the C m hysteresis loops at 
the lower frequencies of 0.12 and 0.05 Hz as with 
the 7=0.43 Hz case. This confirms that the strategy 
of prescribing 360 steps per cycle and between 30 
to 50 sub iterations continues to hold in the low 
frequency range, although the results in Fig. 9 
offer a warning of potential difficulties in 
converging dynamic solutions at lower frequencies. For this particular case, reasonably converged solutions can be 
obtained on the GTM across the desired frequency range with as few as 10800 total solution iterations per pitch 
cycle. Such computations require approximately 275 CPU hours per cycle on an Intel Xeon® 3 GHz cluster for the 
USM3D/SA code. 


Table 3 - GTM cases to evaluate less subiteration at various 


frequencies. USM3D/SA, f=0.43, 0.12, 0.05 Hz, a 0 = 24 deg, 
Aa=+/- 5 deg, 360 steps per cycle. 


Strategy 

Total Iter 

Frequency, Hz 

0.43 

0.12 

0.05 

“Well Converged” 

36000 

360/100 

360/100 

360/100 

“Intermediate” 

18000 

360/50 

360/50 

360/50 

“Most Efficient” 

10800 

360/30 

360/30 

360/30 


Table 2 - GTM cases for determining the 
“most efficient” solution strategy. USM3D/ 
SA,/=0.43 Hz, a o=24 deg, Aa =+/- 5 deg. 


Iterations 
per cycl 

Time steps per 
cycle/ subiteration 

36000 

360/100 

720/50 

1440/25 

18000 

360/50 

720/25 

1440/13 

14400 

360/40 

720/20 

1440/10 

10800 

360/30 

720/15 

1440/8 

7200 

360/20 



3600 

360/10 
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D. Testing Strategy for Advanced Turbulence Model 

The initial stages of the guideline study were 

conducted with the SA turbulence model. Similar Table 4 - GTM cases to evaluate strategy on advanced 
computations are performed with the SST model SST turbulence model at various frequencies, 

to analyze impact of the strategy on a more USM3D/(SST and SA), f=0.43, 0.12, 0.05 Hz, a=24 deg, 

advanced turbulence model. As noted in Table 4, Aa=+/- 5 deg, 360 steps per cycle, 

the computations are performed for a 0 = 24 deg, 

Aa=+/- 5 deg across the full range of pitch 
frequencies 0.43, 0.12, and 0.05 Hz with 360 
steps per cycle using 30 and 50 subiterations. The 
comparisons in Fig. 13 demonstrate that the 
variation from 30 to 50 subiterations does not 
significantly alter the behavior of the C m 
hysteresis for either of the turbulence models. 

Both turbulence models demonstrate the capacity 
to capture the general behavior of C m hysteresis 
for the GTM at this post-stall condition. While some changes of shape in the C m hysteresis are apparent between the 
SST and SA models, the primary difference is an offset in C m levels that is consistent with the offset in static C m in 
Fig. 6 at a=24 deg. These results highlight the strong impact of turbulence model on the post-stall solution that must 
be considered carefully. 

E. Evaluation for Larger Pitch Amplitude 

The devised strategy is next evaluated for a more 

challenging dynamic condition with larger pitch Tab)e g _ GTM cases t0 evaluate strat on la 

amplitude Ac^+I - 10 deg using 360 time steps per pjtch amplitude . USM 3D/(SST and SA), f=0.43, 0.12, 
cycle. The larger of the two prescribed subiterations 0 Q5 Hz a=24 deg? Aa=+/ _ 10 deg 360 steps per cycle> 
in Table 5 is set to 100 to insure a well-converged 
reference condition for this more difficult condition. 

Both the SA and SST turbulence models are also 
retained for the comparisons. The computations 
presented in Fig. 14 demonstrate that the SA and 
SST solutions each capture the general size of the 
experimental C m hysteresis loops across the 
frequency range. Some differences are apparent in 
the details of the shapes. However as discussed 
earlier, the “mean” wind tunnel data is extracted 
from highly unsteady dynamic measurements and 
its qualitative similarity with the URANS CFD result is remarkable. The actual “goodness” or lack thereof of the 
computational results remains an open question to be resolved. 

F. Evaluation for Grid Sensitivities 

Finally, an assessment of grid sensitivities on the periodic oscillating solution is provided using the two finer 
grids having 12- and 24-million cells. The computations were performed with the USM3D/SA model at y=0.43 Hz, 
a 0 = 24 deg, Aa =+/- 5 deg using 360 steps per cycle and 100 subiterations. As observed in Fig. 15, the general shape 
of the C m hysteresis is maintained for the three grids at y=0.43 Hz, but tend toward collapse at the extrema for the 
two lower frequencies of 0.12 and 0.05 Hz. The overall reduced vertical dimensions of the ellipses properly reflect 
the reduced damping force as frequency (and corresponding maximum pitch rates) is reduced. The small differences 
in C m offset are consistent with the offset of static C m in Fig. 5 at a = 24 deg. In comparison to the analysis from the 
preceding two sections, these results demonstrate that the choice of turbulence model has far more impact on 
solution correlation than the grid density. Hence as a general rule, numerical guidelines can be developed on a 
relatively coarse grid for this class of problem, as long as it has sufficient resolution to produce static solutions to a 
reasonably level of accuracy. 
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IX. Preliminary Harmonic Analysis 

Harmonic analysis offers a common framework to compare dynamic models from two independent modeling 
methodologies. Figures 16 compares linear harmonic models estimated from both wind tunnel and CFD time 
histories. In addition, wind tunnel measurements, mean wind tunnel measurements, and USM3D/SA predictions are 
shown. Three different frequencies of oscillation, at 0.43, 0.12, and 0.05 Hz, with nominal oscillations about a 0 = 24 
degrees and with Aa =+/- 5 degrees were analyzed. In this high angle of attack regime, well beyond stall, nonlinear 
and unsteady behaviors are expected, however the relatively small amplitude of oscillation limits the degree of 
nonlinear responses. 

The deviation of measured and computed hysteresis loops from the 1 st harmonic ellipses (dashed lines) in Fig. 16 
reflects the degree of aerodynamic nonlinearity. Consequently, the regular ellipses formed by the linear harmonic 
models represent the best fit to the nonlinear responses shown. Excluding bias errors (that reflect static model 
differences) between the two harmonic models, the differences in loop size at the same oscillation frequency reflect 
the differences in the out-of-phase coefficients. Overall matching of the harmonic model responses from the two 
independent modeling methodologies is encouraging and supports further investigation. 

In-phase and out-of-phase coefficient estimates from this analysis were consistent between the two modeling 
methodologies. Out-of-phase coefficients estimated from wind tunnel and USM3D data both reflect a clear 
dependence on frequency and therefore the presence of unsteady behaviors. Differences in coefficient values, 
averaged over 3 frequencies, were 18% for C Nq and 42% for C mq . In-phase coefficients from both sources were also 
consistent and demonstrated negligible frequency dependence. Differences in these coefficient values averaged 1% 
for C Na and 7% for C Wa .. In-phase coefficient values are comparable with the static data slopes. Observing a larger 
difference for both harmonic models in estimating out-of-phase coefficients compared to in-phase coefficients is 
consistent with general aircraft identification experience in estimating static versus dynamic terms. 

Modeling dynamic systems adds a requirement that sufficient sample rates be used that satisfy the Shannon 
Sampling Theorem. Wind tunnel measurements were made with a sample rate of 250 Hz, an anti-alias filter at 100 
Hz, and a final low-pass digital filter at 4 Hz. In this study, sample rates in Hertz for the CFD data were dependent 
upon the product of oscillation frequency and samples/cycle. Consequently, the lowest sample rates of 18 Hz 
occurred for the lowest frequency cases. For this particular study, aircraft configuration, and conditions considered, 
sample rates of approximately 1 8 Hz may be sufficient for the CFD methodology. 

X. Summary 

A systematic study has been presented using the NASA Generic Transport Model to guide the selection of a 
numerical solution strategy for URANS computation of a subsonic transport configuration undergoing simulated 
forced oscillation about its pitch axis. Forced oscillation is central to the prevalent wind tunnel methodology for 
quantifying aircraft dynamic stability derivatives from force and moment coefficients, which is the ultimate goal for 
the computational simulations. A preliminary linear harmonic analysis was also performed to demonstrate the 
potential of extracting reasonable values for dynamic stability derivatives from computational solutions. 

Several key insights can be summarized from the preceding analysis. The following assumes the dynamic 
simulations are performed with a URANS flow solver and turbulence model using a 2 nd -order time stepping scheme, 
and low-speed flow on a subsonic transport configuration. Among the observations noted from the GTM study are: 

• Convergence of periodic oscillating solutions is a function of the total numbers of solution iterations. 
Figures 9 and 10 confirm that for any angle of attack the dynamic solutions converge to the same 
answer regardless of the path as long as (total solution iterations per cycle) = (number of time steps per 
cycle)* (number of subiterations) remain constant. There are apparently no shortcuts to convergence. 

Hence in future attempts to establish time step and subiteration guidelines, there is no need to run a 
wide range of time steps, as in Table 1, for the purpose of determining solution convergence. One 
carefully selected case should be adequate for determining a solution strategy. 

• Dynamic solutions may be more difficult to converge at the lower oscillation frequencies. Subiterative 
solution convergence between time steps n and n+1 at the lower frequencies of 0.12 and 0.05 Hz 
during the GTM forced pitch oscillation cycle exhibit reluctant numerical convergence and an 
unsteady behavior in Fig. 9. This suggests a stronger potential to encounter convergence difficulties at 
lower frequencies than at higher ones. Hence, one should consider starting a convergence guideline 
study by focusing on a “potentially difficult” case at the lowest frequency of interest so that the 
resulting guidelines should readily apply at the better-behaved higher frequencies. 

9 

American Institute of Aeronautics and Astronautics 



AIAA 2010-4819 
Draft 6/10/10 @9:50am 


• A dynamic convergence study can be performed on a reasonably coarse grid to conserve resources. 
Once a solution strategy has been identified, sensitivities to grid resolution should be investigated 
further with the determined strategy. 

• For forced oscillation problems in regions of post-stall flow, the selection of turbulence model has far 
more impact than grid resolution. Ref. [6] demonstrates how the static force and moment coefficients 
can provide useful guidance for dynamic studies. The single a 0 =2A deg condition was selected for the 
present study from Fig. 5 as representative of a post- stall flight condition. 

• Simulations of forced oscillation in post-stall flows with the URANS formulation yields remarkable 
correlation with the time-averaged “ mean ” of highly unsteady dynamic wind tunnel data. The URANS 
formulation assumes Reynolds averaging of the Navier- Stokes equations that effectively filters out any 
unsteady turbulent eddies from massively separated flows. Yet, the overall size and shape of the 
computed C m hysteresis are very similar to the time-averaged “mean” of highly unsteady forced- 
oscillation wind tunnel data. 

• The conclusions from this analysis of a subsonic transport configuration also apply to a fundamentally 
different configuration and flow field. Each of the preceding observations is fully consistent with those 
for a different class of unmanned air combat vehicle dominated by vortical flow in Ref. [6]. 

• The question of “how accurate do dynamic CFD solutions have to be” remains unanswered. Since the 
final product of CFD are dynamic stability derivatives for flight simulator models, further analysis of 
the URANS solutions using harmonic analysis and system identification techniques is warranted 
before pursuing more costly high-resolution turbulence models, such as Detached Eddy Simulation. A 
preliminary harmonic analysis of the GTM solutions suggests that linear out-of-phase dynamic 
stability derivatives can be extracted from CFD within 18- to 42-percent of experimentally derived 
linear derivatives. 

• Largest differences between wind tunnel and CFD harmonic models were observed for pitching 
moment models at lower frequencies. Lower frequencies present a more difficult estimation problem 
for both wind tunnel and CFD cases since the damping forces are smaller leading to reduced 
signal/noise ratio for identification. Future identification work should include confidence interval 
estimation for both CFD and wind tunnel parameter estimates that will allow better accuracy 
comparisons. 

• Frequency dependence was observed in the CFD out-of-phase coefficients. This result reflects the 
presence of unsteady behaviors and is consistent with that obtained from the wind tunnel estimates. 
Overall consistency between the two methodologies in predicting dynamic response supports further 
investigation into applying CFD methods to difficult stability and control modeling problems. 

• Broader input spectrums for both the CFD and wind tunnel testing methodologies are desirable. 
Additional frequencies of sinusoidal forced oscillations tests would allow estimation of more general 
dynamic models. However, advanced identification inputs such wide-band inputs or arbitrary motion 
would allow more efficient estimation of general mathematical models that predict nonlinear responses 
and unsteady aerodynamic lags. 
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Figure 3. Generic Transport Model (GTM) in the NASA Langley 14- by 22-Foot Low-Speed Tunnel. Left 
static mount. Right - pitch oscillation mount. 
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Figure 4. - Surface triangulation for half-span GTM. Volume contains 6-million tetrahedral cells. 



0 20 40 60 80 0 20 40 60 80 

a, deg. a, deg. 


Figure 5. Effect of GTM grid refinement on static lift and pitching moment coefficient. 14X22 71/^=0.077, 
USM3D/SA M USM3D = 0.2, tf<? cre/ =0.54xl0 6 . 
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Figure 6. Effect of turbulence model on static lift and pitching moment coefficient on 6-million cell grid. 
14X22 M„=0.077, USM3D/SA M USM 3d=0.2, tf<? cre y=0.54xl0 6 . 



Figure 7. Raw and mean dynamic pitching moment coefficients from 14X22 wind tunnel. M x = 0.077, 
Re cre f= 0.54x1 0 6 , a 0 =24 deg, Acx=+/-5 deg. 


USM3D time-accurate flow solution 
with pitch oscillation of 5 deg 



Aa=-5° 


Pitching Moment Hysteresis Subiterative Convergence of 

p. t Pitching Moment at Aa = -3.5 degrees 



Figure 8. Approach for establishing visual correlation between time step and subiterative convergence for 
GTM undergoing pitch oscillation. 
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f=0.43 Hz 


f=0. 12 Hz 


f=0.05 Hz 


<D 



Figure 9. Correlation of number of 2 nd -order time steps with subiterations on convergence at 
downswing Aa=- 3.5 deg for GTM. Constant 36,000 total iterations maintained. USM3D/SA M US m3d= 0.2, 
Re cre f= 0.54x1 0 6 , #0=24 deg, Aa =+/- 5 deg. 
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Figure 10. The effects of combined (time-step)*(subiteration)=36000 iterations on GTM dynamic pitching 
moment coefficients. 14X22 71/^=0.077, USM3D/SA M US m3d = 0.2, Re cre f= 0.54xl0 6 , a 0 = 24 deg, Aa =+/- 5 deg. 
L-R: Frequencies of 0.43, 0.12, 0.05 Hz 
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Figure 11. - Convergence correlations of less subiteration wth GTM dynamic pitching moment coefficient at 
/=0.43 Hz. 14X22 M„,=0.077, USM3D/SA M USM 3d= 0.2, /fc cre/ =0.54x10 6 , a 0 = 24 deg, Aa =+/- 5 deg. L-R: 
Number of steps per cycle a) 360, b) 720, and c) 14400 
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Figure 12.- Effect of less subiteration across GTM dynamic frequency range using 360 time steps per cycle. 
14X22 M< = 0.077, USM3D/SA M US m3d= 0.2, Re cre f= 0.54xl0 6 , a^=24 deg, Aa=+I - 5 deg. L-R: Frequencies a) 
0.43 Hz, b) 0.12 Hz, and c) 0.05 Hz. 
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Figure 13. Effect of less subiteration on GTM with advanced SST turbulence model using 360 time steps 
per cycle. 14X22 M= 0.077, USM3D/(SA and SST) M USM3D = 0.2, /?e„ t ,/=0.54x10 6 , a#=24 deg, Aa=+/-5 deg. 
L-R: Frequencies a) 0.43 Hz, b) 0.12 Hz, and c) 0.05 Hz. 





Figure 14. Effect of subiteration on GTM with larger amplitude of pitch oscillation (a 0 = 24 deg, Aa=+I- 10 
deg) using 360 time steps per cycle. 14X22 71/^=0.077, USM3D/(SA and SST) M US m3d= 0.2, Re cre f= 0.54xl0 6 , 
L-R: Frequencies a) 0.43 Hz, b) 0.12 Hz, and c) 0.05 Hz. 



Figure 15. Effect of grid density on GTM dynamic pitching moment coefficient using 360 time steps per 
cycle, 100 subiterations. 14X22 71/^=0.077, USM3D/SA M US m 3 d = ^-^ Re C r e f= 0.54xl0 6 , a 0 = 24 deg, Aa =+/- 5 deg. 
L-R: 6-million, 12-million, and 24-million cells for half-span. 
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Figure 16. Preliminary harmonic analysis on dynamic pitching moment coefficient. 14X22 71/^=0.077, 
USM3D/SA M US m3d = 0.2, 360 steps per cycle, 100 subiterations, Re cre f= 0.54xl0 6 , 0^=24 deg, Aa =+/- 5 deg. L- 
R: Frequencies a) 0.43 Hz, b) 0.12 Hz, and c) 0.05 Hz. 
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